Work distribution functions for hysteresis loops in a single-spin system 
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We compute the distribution of the work done in driving a single Ising spin with a time-dependent 
magnetic field. Using Glauber dynamics we perform Monte-Carlo simulations to find the work 
distributions at different driving rates. We find that in general the work-distributions are broad 
with a significant probability for processes with negative dissipated work . The special ceises of slow 
and fast driving rates are studied analytically. We verify that various work fluctuation theorems 
corresponding to equilibrium initial states are satisfied while a steady state version is not. 
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I. INTRODUCTION 

Consider a magnetic system in a time-dependent mag- 
netic field. Assume that magnetic field is varied peri- 
odically. Then plotting the magnetization of the system 
against the instantaneous magnetic field we get the well- 
known hysteresis curve. The area enclosed by the hys- 
teresis loop gives the work done on the system by the 
external field and this leads to heating of the magnet. In 
the usual picture that one has of hysteresis one expects 
that the work done is positive. However if the magnetic 
system is small (i.e contains small number of magnetic 
moments) then this is no longer true. For a small magnet 
one finds that the work is a fluctuating quantity and in 
a particular realization of the hysteresis experiment one 
could actually find that the magnet cools and does work 
on the driving force. 

In general, for a small system driven by time- 
dependent forces whose rates are not slow compared to 
relaxation times, one typically finds that various nonequi- 
librium quantities, such as the work done or the heat ex- 
changed, take values from a distribution. Recently there 
has been a lot of interest in the properties of such dis- 
tributions. Part of the reason for the interest is that it 
leads us to examine the question as to how the usual laws 
of thermodynamics, which are true for macroscopic sys- 
tems, need to be modified when we deal with mesoscopic 
systems [1]. 

For instance in our example of the magnet with a small 
number of spins, there is a finite probability that all the 
spins could suddenly spontaneously flip against the direc- 
tion of the fleld by drawing energy from the heat bath. 
Intuitively this gives one the feeling that there has been 
a violation of the second law. In fact historically early 
observers of Brownian motion had the same feeling when 
they saw the "perpetual" motion of the Brownian par- 
ticles [2]. However if one looks at the precise statement 
of the second law one realizes that there is no real vio- 
lation. The second law is a statement on the most prob- 
able behavior while here we are looking at fluctuations 
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about the most probable values. These become extremely 
small for thermodynamic systems. On the other hand, 
for small systems these fluctuations are significant and 
a study of the properties of these fluctuations could pro- 
vide us with a better understanding of the meaning of the 
second law in the present context. This will be necessary 
for an understanding of the behavior of mesoscopic sys- 
tems such as molecular motors, nanomagnets, quantum 
dots etc. which are currently areas of active experimental 
interest. 

Much of the recent interest on these nonequilibrium 
fluctuations has focused on two interesting results on the 
distribution of the fluctuations. These are (1) the Jarzyn- 
ski relation [3, 4] and (2) the fluctuation theorems [5-9]. 
A large number of studies, both theoretical and experi- 
mental [10-12] have looked at the validity of these theo- 
rems in a variety of systems and also their implications. 
At a fundamental level both these theorems give some 
measure of "second law violations" . At a practical level 
the possibility of using these theorems to determine the 
equilibrium free energy proflle of systems using data from 
nonequilibrium experiments and simulations has been ex- 
plored [14]. 

In this paper we will be interested in the fluctuations of 
the area under a hysteresis loop for a small magnet. We 
look at the simplest example, namely a single Ising spin 
in a time-dependent magnetic field and evolving through 
Glauber dynamics. Hysteresis in kinetic Ising systems 
have been studied earlier [15, 16] where the main aim 
was to understand various features such as dependence 
of the average loop area on sweeping rates and ampli- 
tudes, system size effects and dynamical phase transi- 
tions . The area distribution was also studied in [16] but 
the emphasis was on different aspects and so is quite in- 
complete from our present viewpoint. A two state model 
with Markovian dynamics was earlier studied by Ritort 
et al [17] to analyze experiments on stretching of sin- 
gle molecules. Systems with more than two states have 
also been studied [18] in the context of single-molecule 
experiments. However the detailed form of the work- 
distributions have not been investigated and that is the 
main aim of this paper. These distributions are of inter- 
est since there are only a few examples where the explicit 
forms of the distributions have actually been worked out 
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[17 20]. Most of the experiments, for example those in 
the RNA stretching experiments of Liphardt et al [10] 
or the more recent experiment of Douarche et al [13] on 
torsionally driven mirrors, are in regimes where the work- 
distributions are Gaussian. 

We perform Monte-Carlo simulations to obtain the dis- 
tributions for different driving rates. We consider differ- 
ent driving protocols and look at the two cases corre- 
sponding to the transient and the steady state fluctu- 
ation theorems. It is shown that the limiting cases of 
slow and fast driving rates can be solved analytically. 
We also point out that the problem of computing work- 
distributions is similar to that of computing residence- 
time distributions. 

The organization of the paper is the following: In sec- 
tion (II) we define the dynamics and derive the Fokker- 
Planck equation satisfied by the distribution function. In 
section (III) we give the results from Monte-Carlo simu- 
lations and in section (IV) we treat the special cases of 
slow and fast driving rates. Finally in section (V) we give 
a discussion of the results. 



II. DEFINITION OF MODEL AND DYNAMICS 

Consider a single spin, with magnetic moment fi, in a 
time-dependent magnetic field h{t). The Hamiltonian is 
given by 



H = —fiha (7 = ±1 



(1) 



We assume that the time-evolution of the spin is given by 
the Glauber dynamics. Let us first consider a discretized 
version of the dynamics. Let the value of the magnetic 
field at the end of the (n — l)th time step be and 
let the value of the spin be (Tn-i. The discrete dynamics 
consists of two distinct process during the nth time step: 

1. The field is changed from to /i„ = /i„_i-|-A/i„. 
During this step an amount of work AW = — /i(T„_iA/i„ 
is done on the system. 

2. The spin flips with probability ^(6"'^'*'''''^"- V-^) 
where Z = e'^^'*" -|- e'^^'^" is the equilibrium partition 
function at the instantaneous field value. The factor p is 
a parameter that is required when we take the continuum 
time limit and whose value will set equilibration times. 
At the end of this step the spin is in the state C7„. During 
this step the system takes in an amount of heat AQ = 
—fihnicrn — cr„_i) from the reservoir. 

Given the microscopic dynamics we can derive time- 
evolution equations for various probability distributions. 
These are standard results but we reproduce them here 
for completeness. 

Fokker- Planck equation for spin distribution: First 
let us consider the spin configuration probability Pn{o') 
which gives the probability that at time n the spin is in 
the state a. We write the field in the form ft,„ = ft-o/n 
where /„ is dimensionless and let us define e = P/iho. 



Then we get the following evolution equation: 



p- 



To go to the continuum-time limit we take the limits 
p^O, At^Q, with p/At r and /„ f{t), Pn{cr) 
P{(T, t). Using the dimensionless time t = rtwe then get: 



dP 
a7 



= —TP where 



(2) 



^ - ' P(i,r) j' ^ 



Z 



The magnetization m(r) = 2P(|, r) — 1 thus satisfies the 
equation 



dm{T) 
dr 

whose solution is 



-m(T) + tan/i[e/(r)] 



(3) 



m(T) = e 



-^m{0)+ f dT'e"(^^^'Uanh[e/(T')]. (4) 

^0 



Fokker- Planck equaMon for Work distribution: The to- 
tal work done at the end of the nth time step is given by: 



-j»y~^o"„_iA/i„ 



(5) 



To write evolution equations for the work- distribution it 

is necessary to first define Qn{W, cr), the joint probability 
that at the end of the nth step the spin is in the state a 
and the total work done on it is W. Then Qn{W,a) will 
satisfy the following recursions: 

Qn+l (W,t) = {l-p^-^))Qn{W + i,Ahn,^) 



ri+1 



+ )Qn{W - llAhnA) 

Qn+i{w,i) = p(^^^)Q„(I'F + mA/i„,T) 

+ {l-p{^^))Qn{W-flAhn,l). 

Taking limits At ^ 0,p ^ with p/ (At) -^r,hn^ h{t) 
and Ahn/{At) h and using dimensionless variables r 
and w = (3W = — e /J"" dTm{T)df / dr we finally get 

— = -T.Q + e—a^.— where (6) 
OT dr aw 



_ . Q{w,],t) 
^ ~ ' Q{w,i,r) j ' 



1 
-1 



From Eq. 6 we get the following equation for Q{w,t) 
Q{w,tT) + Q{w,l,T): 



+ (1 - = e/tanh(e/)|^ + (e/)^^'^ 



f'dr 



dw 



(7) 
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We have not been able to solve these equations ana- 
lytically except in the limiting cases where the rate of 
change of the magnetic field is very slow or very fast. In 
the next two sections we will first present results from 
Monte-Carlo simulations which give accurate results for 
any rates and then discuss the special cases. 

III. RESULTS FROM MONTE-CARLO 
SIMULATIONS 



In the following sections we will verify that Crook's FT 
is always satisfied but that TFT does not hold whenever 
the distributions are non-Gaussian processes and the pro- 
cess does not have time-reversal symmetry. We will also 
test the validity of the steady state fluctuation theorem 
(SSFT). This has the same form as Eq. 9 with the differ- 
ence that the initial state is chosen from a nonequilibrium 
steady state distribution instead of an equilibrium distri- 
bution. 



We have studied three different driving processes: 

(A) The system is initially in equilibrium at zero field 
and the field is then increased linearly from to e 
in time Tm- 

(B) The system is initially in equilibrium and the field, 
which is taken to be piecewise linear, is changed over one 
cycle. 

(C) The system is run through many cycles till it 
reaches a nonequilibrium steady state. We measure work 
fluctuations in this steady state. 

In cases (A) and (B) we will be interested in testing the 
transient fluctuation theorem (TFT) while in case (C) we 
will look at the steady state fluctuation theorem (SSFT). 
Let us briefly recall the statements of these theorems for 
work distributions in systems with Markovian dynamics. 

Crook's Fluctuation Theorem: In this case the system 
is initially in thermal equilibrium and then an external 
parameter [e.g magnetic field h{t)\ is changed from an 
initial value hi, at time t = 0, to a final value /i/ in 
a finite time tm- Suppose the work done on the sys- 
tem during this process is W and the change in equi- 
librium free energy is AF. Prom thermodynamics we 
expect that for a quasistatic process W = AF. Hence 
for our nonequilibrium process it is natural to define the 
dissipated work as Wd = W — AF which has a distribu- 
tion Q{Wd)- Now consider a time-reversed path for the 
external field /i./j(t) = h{tm — t) for which the work distri- 
bution is QRiWd). The fluctuation theorem of Crook's 
then states: 



Q{Wd) 



(8) 



For Gaussian processes it can be shown that QR{Wd) = 
Q{Wd) [8] and hence we get the usual form of the tran- 
sient fluctuation theorem (TFT) 



QjWd) 
Q{-Wd) 



(9) 



Another situation where TFT is satisfied is the case 
where the field is kept constant or if the process is time- 
reversal symmetric. Finally we note that the Jarzynski 
relation 



^e-0w,^ = 1 



(10) 



follows immediately from Crook's theorem Eq. 8 and so 
will be satisfied in all cases where Crook's holds. 



A. Field increased linearly from to e 

In this case /(r) = r/Tm and we have chosen e = 0.5. 
We note that with a static field, the equilibration time 
is given by = l/r or = 1. The rate of change of 
magnetic field is 1 /rm and slow and fast rates correspond 
to large and small values for r,„ respectively. In Fig. 1 
we plot the work distributions for various values of Tm- 
We have plotted the distribution of the dissipated work 
Wd ^ w - 13 AF (Here /3AF = - In cosh e). In Fig. 2 
we plot the average magnetization as a function of field, 
again for different rates. Some interesting features of the 
work-distributions are: 

(i) The distributions are in general broad. This is true 
even at the slowest driving rates where the average mag- 
netization (Fig. 2) itself is close to the equilibrium pre- 
diction. Note that the allowed range of values of Wd is 
[-e - (5AF, e - PAF] w [-0.38, 0.62]. Also we see that 
the probability of negative dissipated work is significant. 

(ii) For slow rates the distributions are Gaussian and 
this can be understood in the following way. Imagine 
dividing the time range into small intervals. Because 
the rate is slow, there are a large number of spin flips 
within each suc;li interval, and so the average magnetiza- 
tion from one interval to the next can be expected to be 
uncorrelated. Since the work is a weighted sum of the 
magnetization over all the time intervals we can expect 
it to be a Gaussian. 

(iii) For fast rates we get 5— function peaks at w = ±e. 
This again is easy to understand since the spin doesn't 
have time to react and stays in its initial state. In section 
(IV) we will work out analytic expressions for the the 
work distributions by considering probabilities of 0— spin 
flip and 1— spin flip processes. 

For slow rates we have verifled (see Fig. 3) that the 
fluctuation theorem is satisfied. For faster rates we see 
that the probability of negative work processes is higher 
than what is predicted by the TFT. 



B. Field is taken around a cycle 

As shown in Fig. 4 we consider two different cyclic 
forms for /(r). One is a symmetric cycle and the 
other a asymmetric one. For these two cases the work- 
distributions are plotted in Fig. 5 and Fig. 6 respectively. 
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FIG. 1: Distributions of the work done in driving a magnet 
at different rates wlien magnetic field is changed linearly. 
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FIG. 3: Plot shows that the fluctuation theorem is valid for 
slow processes with Gaussian work-distributions. Inset shows 
that for a fast rate (rm = 5.0) the probability of negative 
work is much larger than that predicted by the FT. 
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FIG. 2: Average magnetization m(r) for different rates, when 
magnetic field is changed linearly. 

For the symmetric cycle we plot the average magnetiza- 
tion as a function of the field in Fig. 7. This gives the 
familiar hysteresis curves. 

As before we again find that the work-distributions 
are broad. For slow rates we get Gaussian distributions 
while for fast rates we get a (5— function peak at the origin 
which correspond to a 0— spin flip process. The slow and 
fast cases are treated analytically in section (III). 

As expected we can verify the transient fluctuation the- 
orem for both the symmetric and asymmetric processes. 
That TFT should be satisfied follows from Crooks FT 
and noting that the time reversed process has the same 
distribution as the forward process because of the addi- 
tional ft, — > — ft, symmetry that we have in this case. 

We have also studied an asymmetric half- cycle for 




FIG. 4: Magnetic field changed over a cycle. In the symmetric 
case the total cycle time is 4rm while in the asymmetric case 

it is 2{Tm +Tn). 



which Qniwd) Q{wd)- Consequently we find that the 
usual TFT is not satisfied while the more general form 
of TFT of Crooks holds. We show this in Fig. 8 where 
we have plotted Q{wd), e^'^Q{—Wd) and e'^'^QR{—Wd)- 



C. Properties in the nonequilibrium steady state 

We now look at the case when the spin is driven by the 
oscillating field into a nonequilibrium steady state and 
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FIG. 6: Plot of work-distributions obtained for different asym- 
metric cycles of the magnetic field. 

we measure fluctuations in this steady state. In this case 
the work distributions (over a cycle) have the same forms 
as in the transient case (Fig. 9). The joint distribution 
function Q{w,a,T) satisfies the same equation Eq. 6 but 
now the initial conditions are different. In Fig. 10 we plot 
the steady state hysteresis curves. Note that unlike the 
transient case the hysteresis curves are now closed loops. 

Finally we test the validity of the steady state fluctua- 
tion theorem (SSFT). This theorem has been proved for 
dynamical systems evolving through deterministic equa- 
tions but there exists no proof that a similar result holds 
for stochastic dynamics. From Fig. 11 it is clear that 
SSFT does not hold. 
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FIG. 8: Plot showing the work-distribution for an asymmet- 
ric half-cycle and the validity of Crook's fluctuation theorem. 
Note that the probability of negative work processes is much 
higher than that predicted by usual TFT. 



IV. ANALYTIC RESULTS FOR SLOW AND 
FAST RATES 

A. Field increased linearly from to e 

(i) Slow case : >> 1 

As argued in the previous section we expect that the 
work-distributions to be Gaussians which will be of the 
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FIG. 9: Work-distributions in the nonequilibrium steady pj^ Violation of the fluctuation theorem for the steady 

state work-distribution corresponding to Tm = 0.5. Inset 
shows the same plots for the transient case where FT is clearly 
satisfied. 
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FIG. 10: Hysteresis curves in the nonequilibrium steady state. 



general form 



Qiw) 



1 



(11) 



Since the distribution satisfies Jarzynski's equality, it fol- 
lows at once that the mean and variance are related by 



(12) 



Hence we just need to find the mean work done. The 
mean work done is given by (w) = —{e/Tm) Jq"^ dTm{T). 
In the strict adiabatic limit r„i — s- oo we have 
nT-adiT) = tanh{eT/Tm) and the mean work done {w) = 



— log(cosh(e)) — PAF. For large we try the pertur- 
bative solution 



1 



TO(r) = madir) H 5(t) 



(13) 



Substituting in Eq. 3 we get an equation for ^(t) whose 
solution gives 



For the work done we then get 



(3AF 



-tanh{e) 



(14) 



(15) 



In Fig. 12 we compare the simulations for slow rates with 
the analytic results. 

(u) Fast case : Tm « 1- 

If we change the field very fast then the spin is not 
able to respond and so there are few spin flips during the 
entire process. At the lowest order there is no flip and 
this gives rise to the 5— functions peaks at ±e seen in the 
distribution. We now calculate the work-distribution by 
looking at contributions from 0— spin flip and 1— spin flip 
processes. Let S('\,to,t) be the probability that, given 
that the spin is t at time tq, it remains in the same state 
till time r. It is easy to see that 5(t,To,r) satisfies the 
equation 



-5(t,ro,r) 



Dt Z 
Solving we get, for the linear case /(r) ~ t/t„ 



(16) 



5(T,0,r) 



Jo 



■^I^'^"' =e^(cosh( — ))* (17) 
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FIG. 12: Comparison of work-distributions for slow rates ob- FIG. 13: Comparison of work-distribution for a fast rate ob- 
tained from simulations and from the analytic form. Solid Gained from simulation and from the analytic form, 
lines show the analytic results. 



Putting T = T„i corresponds to the process for which the 
work done is w = — e. Hence, since the probabihty of the 
spin being initially in "f state is 1/2, we get 

1 

Prob{w ^ -e) = [cosh(e)]— (18) 

Proceeding in a similar fashion by starting with a | spin 
we get 

Prob{w = = ^ e~?^(cosh(e))~^ (19) 

Next let us consider 1— spin flip processes which (for fast 
rates) are the major contributors to the part of the dis- 
tribution between the two peaks. Let Si{'\,T)dT be the 
probability that the spin starts in the "f state, flips once 
between times t to t + dr, and stays J, till time Tm- This 
is given by 

^i(T, T)dT = Sit 0, T)(e-^^/^'"/Z)dT5a, T, r„) (20) 
The work done during such a process is given by 

w = -—{2t-t^) (21) 

Similarly the case where the spin starts from a J, state 
gives 

Sl{l,T)dT = 5(i,0,T)(e^^/^"VZ)dT5(T,T,T„0 (22) 
and the work done in this case is 

w^ — {2t- t™) (23) 



Adding this two contributions and plugging in the form 
of S{a, To, t) obtained earlier we get the following contri- 
bution to the work-distribution: 

TVi -i^rr,-^) e^(cosh(£)^) et(cosh(£)^) 
^'^""^"86^ ^ cosh(^) ^ cosh(^) 

The full distribution is given by 

Q{w) = ^e~~2^[cosh{e)]^ 6{w + e) 

+ ie"* [cosh(e)]-*5(u; - e) + Qi{w) (24) 

for —e<w<e, and zero elsewhere. In Fig. 13 we 
show a comparison of this analytic form with simulation 
results for Tm — 0.01. The strengths of the 5— functions 
at u; = ±e are accurately given by Eqs. 18, 19. 



B. Field is taken around a cycle 

(i) Slow case : r„i >> 1 

We again expect a Gaussian distribution and since 
AF = for a cyclic process, hence the mean and vari- 
ance of the distribution are related by crj — 2{w) (see 
Fig. 14). As before we compute the mean work to order 
1/t,„ and find 

4e 

(w) = — tanh{e) (25) 

(ii) Fast case : r„i << 1. 

In this case the work distribution gives a (5— function 
peak at the origin for 0— spin flip processes. To find the 
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This has the solution 
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FIG. 14: Comparison of work-distributions for slow rates ob- 
tained from simulations of the cyclic case and from the ana- 
lytic form. Solid lines show the analytic results. 
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FIG. 15: Comparison of work-distribution, for the symmetric 
cycle with a fast rate, obtained from simulations and from the 
analytic results. 



S'(t,0, r„) = e' 



(26) 



Adding up an equal contribution from S'(i,0,Tm), and 
since both initial conditions occur with probability half, 
we finally get 



Prob{w — 0) = e 



-2t„ 



(27) 



Next we look at the contribution of 1— spin flip processes. 
Let the spin flip occur between times r and r + dr. It is 
convenient to divide the total time 4r,„ into four equal 
intervals, the dependence of w on r being different in each 
of the intervals. Thus if we start with the spin initially 
in an "f state then we have 



2er 



w 



<T <T„ 



2e 



(r - 2r„) Trn <T < 2t„ 



2e , 

= (2Tm - t) 2T,n < T < 3Tm 

2e , 

= (r - 4Tm) 3Tm < T < 4Tm- 



The probabilities of each of these processes is again given 
by: 

Si{],T)dT = 5(T,0,r)(e-^/M/Z)dT5a,T,4T™) (28) 

Using the relations between w and r and summing up 
the four different possibilities we finally get (for initial 
spin state t) 



-e-2-'"et[(cosh(^)) 



2 

+ (cosh(|))-^-i(cosh(e))- 



(29) 



for — 1 < w < 1 and zero elsewhere. Note that the al- 
lowed range of w is [—2, 2] but single spin-fiip processes 
only contribute to work in the range [—1, 1]. Similarly if 
we start with spin state J, we get 



probability of this, we solve Eq. 16 with /(t) for the cycle 
given by 



,W . 



e2[(cosh(-)) 

Ot Z 

+ (cosh(^))^-i(cosh(6))- 



(30) 



for — 1 < w < 1. The full work-distribution (contribution 
from 1— spin flip processes) is thus: 



/(t) = — T for < T < T„ 



/(r) = — (2r„ - t) for r„i < t < 3r„ 
/(r) = — (r-4r™) for 3t™ < t < 4r„ 



Q{w) 



-2t„ 



5{w) + Q\{w) + Q{{w) 



(31) 



In Fig. 15 we compare the analytic and simulation results. 
The strength of the (5— function at it; = is accurately 
given by Eq. 27. 
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V. DISCUSSION 

We have computed probability distributions of the 

work done when a single spin, with Markovian dynam- 
ics, is driven by a time dependent magnetic field. We 
find that work fluctuations are quite large (even for slow 
driving rates) and there is significant probability for pro- 
cesses with negative dissipated work. For slow driving the 
number of spin flips during the entire process is very large 
and the total work is effectively a sum of random vari- 
ables. Hence the distributions are Gaussian with widths 
proportional to the driving rate. On the other hand for 
very fast driving the probability of flipping is low and we 
can compute the work-distributions perturbatively from 
probabilities of zero-flip, one-flip, etc. processes. 

While the two special cases of slow and fast rates can be 
solved, it looks difficult to obtain a general solution valid 
for all rates even in this single particle problem. We note 
that the problem of calculating the work-distribution is 
similar to that of calculating residence-time distributions 
in stochastic processes [21-23]. In fact for the case in 
section(lllA) the work done is proportional to the av- 
erage magnetization which is easily related to the resi- 
dence time (time spin spends in t state). For stationary 



stochastic processes, such as the random walk, the resi- 
dence time distribution can be obtained exactly. However 
for non-stationary processes this becomes difficult and no 
exact solutions are available [23] . fn our spin-problem too 
it appears that the non-stationarity of the process makes 
an exact solution difficult. 

For a system with A'^ spins the total work done on the 
system is simply a sum of the work done on each of the 
spins. For the case where the spins are non-interacting we 
thus get a sum of N independent random variables. For 
large N the distribution will be a Gaussian with a mean 
that scales with N and variance as N^^"^. For interacting 
spins the properties of the work-distribution is an open 
problem. Especially of interest is the question as to what 
happens as we cross the transition temperature. Finally 
we note that the large fluctuations in the area under a 
hysteresis curve should be experimentally observable in 
nano-scale magnets. 
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